Introduction to the Diffusion Monte Carlo Method 
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I. INTRODUCTION 

The Schrodinger equation provides the accepted de- 
scription for microscopic phenomena at non-relativistic 
energies. Many molecular and solid state systems 
are governed by this equation. Unfortunately, the 
Schrodinger equation can be solved analytically only in a 
few highly idealized cases; for most realistic systems one 
needs to resort to numerical descriptions. In this paper 
we want to introduce the reader to a relatively recent nu- 
merical method of solving the Schrodinger equation, the 
Diffusion Monte Carlo (DMC) method. This method is 
suitable to describe the ground state of many quantum 
systems. 

The solution of the time-dependent Schrodinger equa- 
tion can be written as a linear superposition of stationary 
states in which the time dependence is given by a phase 
factor exTp(—iE n t/h), where E n is the n-th energy level 
of the quantum system in question. The energy scale 
can be chosen such that all energies are positive. In the 
DMC method one actually considers the solution of the 
Schrodinger equation assuming imaginary time r, i.e., 
after replacing the time t by —it. The solution is then 
given by a sum of transients of the form exp (—E n T/h), 
n = 0, 1, . . .. The DMC method is based upon the ob- 
servation that, as a quantum system evolves in imag- 
inary time, the longest lasting transient corresponds to 
the ground state with energy Eq < E n , n = 1, 2, Fol- 
lowing the evolution of the wave function in imaginary 
time long enough one can determine both the ground 
state energy E and the ground state wave function 4> 
of a quantum system, regardless of the initial state in 
which the system had been prepared. The DMC method 
provides a practical way of evolving in imaginary time 
the wave function of a quantum system and obtaining, 
ultimately, the ground state energy and wave function. 

The DMC method can be formulated in two different 
ways. The first one is based on the similarity between 
the imaginary time Schrodinger equation and a general- 
ized diffusion equation. The kinetic (potential) energy 
term of the Schrodinger equation corresponds to the dif- 
fusion (source/sink or reaction) term in the generalized 



diffusion equation. The diffusion-reaction equation aris- 
ing can be solved by employing stochastic^-palculus as it 
was first suggested by Fermi around 1945Era. Indeed, the 
imaginary time Schrodinger equation can be solved by 
simulating random walks of particles which are subject to 
birth/death processes imposed by the source/sink term. 
The probability distribution of the random walks is iden- 
tical to the wave function. This is possible only for wave 
functions which are positive everywhere, a feature, which 
limits the range of applicability of the DMC method. 
Such a formulation of the, DMC method was given for 
the first time by Andersoro who used this method to cal- 
culate the ground state energy of small molecules such as 

A second formulation of the DMC method arises 
from the Feynman path integral solution of the time- 
dependent Schrodinger equation. By means of path in- 
tegrals the wave function can be expressed as a multi- 
dimensional integral which can be evaluated by employ- 
ing the Monte Carlo method. Algorithms to solve the 
diffusion-reaction equation obeyed by the wave function 
and algorithms to evaluate the path integral representa- 
tion of the wave function yield essentially one and the 
same formulation of the DMC method. Which one of 
the two formulations of the DMC method one adopts 
depends on one's expertise: a formulation of the DMC 
method based on the diffusion-reaction equation requires 
basic knowledge of the theory of stochastic processes; a 
path integral formulation obviously requires familiarity 
with the corresponding formulation of quantum mechan- 
ics. 

The purpose of this article is to provide a self-contained 
and tutorial presentation of the path-integral formula- 
tion of the DMC method. We also present a numer- 
ical algorithm and a computer program based on the 
DMC method and we apply this program to calculate the 
ground state energy and wave function for some sample 
quantum systems. 

The article is organized as follows: The formulation of 
the DMC method is presented in Sec. |l| 



III a 



In Sec 

numerical algorithm for the DMC method is constructed 
The results of the DMC simulation for some simple quan- 



turn mechanical systems are presented in Sec. [V. Fi- 
nally, Sec. |V| provides suggestions for further numerical 
experiments and guides the reader to the literature on 
the DMC method. 



II. THEORY 

The theoretical formulation of the DMC method, pre- 
sented below, follows three steps. These steps will be 
outlined first, to provide the reader with an overview. 
First Step: Imaginary Time Schrodinger Equa- 
tion. In this step, the solution of the time-dependent 
Schrodinger equation of a quantum system is expressed 
as a formal series expansion in terms of the eigenfunctions 
of the Hamiltonian. One then performs a transformation 
from real time t to imaginary time r, replacing t — > —it. 
The solution of the obtained imaginary time Schrodinger 
equation becomes a series of transients which decay ex- 
ponentially as r — ► oo. The longest lasting transient 
corresponds to the ground state (i.e., the state with the 
lowest possible energy) of the system. 
Second Step: Path Integral Formulation and 
Monte Carlo Integration. In this step, the imagi- 
nary time Schrodinger equation is investigated by means 
of the path integral method. By using path integrals 
the solution of this equation can be reduced to quadra- 
ture, provided that an initial state wave function is given. 
Standard Monte Carlo methodsa permit one to evaluate 
numerically the path integral to any desired accuracy, 
assuming that the initial state wave function and, there- 
fore, the ground state wave function as well, is positive 
definite. In this case the wave function itself can be inter- 
preted as a probability density and the "classical" Monte 
Carlo method can be applied. According to the gen- 
eral principles of quantum mechanics, only the square of 
the absolute value of the wave function has the meaning 
of a probability density; the fact that the ground state 
wave function has to be a positive definite real quan- 
tity imposes severe limitations on the applicability of 
the Monte Carlo technique for solving the Schrodinger 
equation. An efficient implementation of the standard 
Monte Carlo algorithm for calculating the wave function 
as a large multi-dimensional integral is realized through 
an alternation of diffusive displacements and of so-called 
birth-death processes applied to a set of imaginary par- 
ticles, termed "replicas" , distributed in the configuration 
space of the system. The spatial distribution of these 
replicas converges to a probability density which repre- 
sents the ground state wave function. The diffusive dis- 
placements and birth-death processes can be simulated 
on a computer using random number generators. 
Third Step: Continuous Estimate of the Ground 
State Energy and Sampling of the Ground State 
Wave Function. In this step, the ground state energy 
and the ground state wave function are actually deter- 
mined. As mentioned above, the Monte Carlo method 



samples the wave function after each time step. The 
spatial coordinate distribution of the replicas involved 
in the combined diffusion and birth-death processes, af- 
ter each (finite) time step, provides an approximation 
to the wave function of the system at that given time. 
The wave function converges in (imaginary) time towards 
the (time-independent) ground state wave function, if 
and only if the origin of the energy scale is equal to the 
ground state energy. Since the ground state energy is ini- 
tially unknown, one starts with a reasonable guess and, 
after each time step in which a diffusive displacement 
and birth-death process is applied to all particles once, 
one improves the estimate of the ground state energy. 
Ultimately, this estimate converges towards the desired 
ground state energy and the distribution of particles con- 
verges to the ground state wave function. 

In the following, we shall provide a detailed account of 
the above steps. 



A. Imaginary Time Schrodinger Equation 

For simplicity, let us consider a single particle of mass 
m which moves along the x-axis in a potential V(x). Its 
wave function ^(x, £)_is governed by the time-dependent 
Schrodinger equationd 



<9* 
where the Hamiltonian has the form 



H 



2m dx 2 



V{x) 



(2.1) 



(2.2) 



Assuming that the potential for x — > ±oo becomes infi- 
nite, i.e., the particle motion is c onfi ned to a finite spatial 
domain, the formal solution of (2.1) can be written as a 
series expansion in terms of the eigenfunctions of H 



\J>(x,£) = yiCn<f> n (x) e R 



E~t 



(2.3) 



The eigenfunctions <f> n (x), which are square-integrable 
in the present case, and the eigenvalues E n are obtained 
from the time-independent Schrodinger equation 



H<j> n (x) = E n cf> n (x) 



(2.4) 



subject to the boundary conditions lim^^ioo 0„ (a;) = 0. 
We label the energy eigenstates byn = 0,1,2,... and 
order the energies 



Eq < Ei < E2 5* 



(2.5) 



The eigenfunctions <p n (x) are assumed to be orthonormal 
and real, i.e., 



dxcj) n (x)(j> m (x) 



(2.6) 



The expansion coefficients c n in (2.3) are then 

dx(/) n (x)^(x,0) , n = 0,1,2,... , (2.7) 



i.e., they describe the overlap of the initial state &( x, 0) , 
else assumed real, with the eigenfunctions 4> n (x) in (2.4). 
Shift of Energy Scale. We perform now a trivial, 
but methodologically crucial shift of the energy scale in- 
troducing the replacements V(x) — > V(x) — Er and 
E n — v EVi — •E'iJ- This leads to the Schrodinger equation 



and to the expansion 



*(#,£) = y^c K K (x) e" 



ia< 



(2.8) 



(2.9) 



n=0 



Wick Rotation of Time. Now let us perform a 
transformation from real time to imaginary time (also 
known as Wick rotation) by introducing the new vari- 
able t — it. The Schrodinger equation (|2.8|) becomes 



d^ 



n 2 s 2 * 



dr 2m dx 2 

and the expansion ([2.9]) reads 



- [V{x) - E R ] 4< 



*(x,t) 



n=0 



C„ 0n(#) e ~ 



(2.10) 



(2.11) 



No ting the energy ordering ( |2.5| ), one can infer from 
( [2.11 ) the following asymptotic behavior for r — > oo 



(i) if Er > £?o, limT-^oo ^(x, r) = oo, the wave func- 
tion diverges exponentially fast; 

(ii) if En < Eq, lim-r^oo ^&(x, r) = 0, the wave func- 
tion vanishes exponentially fast; 

(hi) if En = Eq, lim T ^oo ^{x, t) = Cq4>q(x), the wave 
function conve rges , up to a constant factor cq de- 
fined through (2.7), to the ground state wave func- 
tion. 

This behavior provides the basis of the DMC method. 
For Er = E , the function &(x, t) converges to the 
ground state wave function <fio(x) regardless of the choice 
of the initial wave function ^(a;,0), as long as there is 
a numerically significant overlap between ^(x,0) and 
cj>o(x), i.e., as long as Co is not too small. The ground 
state wave function, for a single particle, has no nodes (in 
case of many fermion systems this might not be true) and 
one can always fulfill the requirement of non-vanishing 
Co by choosing a positive definite initial wave function 
centered in a region of space where <po(x) is sufficiently 
large. 

We now seek a practical way to integrate equation 
(2.10) for an arbitrary reference energy Er and initial 
wave function ^(x, 0). We shall accomplish this by using 
the path integral formalism. 



B. Path Integral Formalism 



The solution of the imaginary time Schrodinger equation 
fl2.10| ) can be written 



}&(x,t) 



dx K(x,T\x ,0)^{xo,0) , (2.12) 



where the propagator K(x, t\xh, 0) is expressed in terms 
of the well-known path integrals, modified by the replace- 
ment t = — IT 



K(x,t\xq,0) = lim / dx\ 

N—>CG 



N 



(2.13) 



At tt-^, \ m 2 „, 1 t-, 

{ ' xp ^ — j~ 1^ [~2Kt^ ^ 3 ~ X} ~ 1 ' + ^ Xj > ~ R 
3=1 



Here At = t/N is a small time step. One sets xn = x. 
The wave function ^(x,t) can be written in the form 



^(x, t) = lim 

N->o> 



poo (n-\ \ N 

/ U dx n n^^™) 

J -°° \3=0 J n=l 



X P(x n ,X n -i)$(x ,0) , 

where we have defined 



(2.14) 



/ m \ 2 
P (**,**-i) = {^^) exp 



m(x n - x n -\Y 



2TiAt 



and 



W (x„) = cxp 



[V (x n ) - E R ] At 



(2.15) 



(2.16) 



The fun ction P (x n , X n -\) is related to the kinetic energy 
term in (2.2). This function can be thought of as a Gaus- 
sian probability density for the random variable x n with 
mean equal to x n — 1 and variance 



a = ^HAt/j 



(2.17) 



The so-called weight function W (x n ) depends on both 
the potential energy in ( |2.2| ) and the reference energy 
Er. The main difference between the functions P and 
W is that the former can be interpreted as a probability 
density since 



dyP(x,y) = 1 



(2.18) 



while the latter can not. 



The path integral (2.14) can be evaluated analyti- 
cally|-pnly for particular forms of the potential energy 
V{xp. Fortu nately , by choosing N sufficiently large, one 
can evaluate ( 2.14 ) numerically to any desired accuracy. 
However, since a suitable N is necessarily a large number, 
the standard algorithms of numerical integration cannot 
be employed directly:, instead, one uses the so-called 
Monte Carlo methodQ. According to this method any 
(convergent) TV-dimensional integral of the form 



,00 (n-i \ 

1 = 1Z[ dx A f( x O,---,XN-l)'P(xo, 

■> ~°° \j=0 J 

where V is a probability density, i.e., 

V (xq, . . . ,xn-i) > 0, and 

/ II dx i ^ ( x °> ' ' ■ ' x N-i) = 1 i 

J -°° \j=o J 



,X N -i) 

(2.19) 



(2.20) 



can be approximated by the expression 
— V f(x® x W 



2 



(2.21) 



i— 1 



Here the notation xW g 7? means that the numbers x^- , 
i = l,2,... , A/"; J = 0, . . . , N — 1, are selected randomly 
with the probability density V . The larger Af, the bet- 
ter is the approximation-^ « X. In fact, according to 
the central limit theoremu, the values of X obtained as 
a result of different simulations are distributed normally, 
i.e., according to a Gaussian distribution, around the ex- 
act value /, the standard deviation being proportional to 

l/VM. 



In order to evaluate \&(x,t) in ( 2.14 ), for given x, r 
and N, one defines 



JV 



V(x , 



. ,XJV- 



-l) = IJ P(Xn,X n -l) , 



and 



N 



/(so,- 



,a;jv. 



.1) = *(x ,o) n^w 



(2.22) 



(2.23) 



n=l 



such that one can apply ( 2.19 ) and ( 2.21 ). For this pur- 
pose we note that due to 



dyP{x,y)P(y,z) = P(z,z) 



(2.24) 



and ( 2.1S| ) the prob abilit y distribution ( 2.22) do es indeed 
obey the property ( 2.20 ). The expression ( |2.21 ) can now 
be in voked to evaluate <3/(x, r) by means the path inte- 
gral (2.14). This requires the generation of sets of co- 



ordinate vectors x^ S V, i^ 



1,2, . . . ,A/\ where x 



K- 






,«$! for 



In order to obtain vectors x*- 1 ' which sample the prob- 
ability density V one proceeds as follows: In a first step 
one generates, for a fixed x = x]y\ a Gaussian random 
number x ] y_ 1 with mean value x^f (i.e., a Gaussian ran- 
dom number distributed about x^ ) and variance u, ac- 
cording to the probability density P \X^ ,x^_ 1 \ given 
by (|2.15). In a second step, a Gaussian random num- 



ber ij_ ,, with mean xjj_ l and variance a, is generated 
according to P ( x]^_ l7 Xjy_ 2 ) ■ The steps are continued 

to produce random numbers x« until one reaches Xq . 

Two consecutive random numbers x„ and x„_i are re- 
lated through the equation 



,« _ ~W 



crp 



(i) 



(2.25) 



where u is given by (2.17) and p n is a Gaussian ran- 



» 



dom number with zero mean and a variance equal to 

(i) 

one. The p n 's can be generated numerically by means 
of algorithms referred t o as rando m nu mber generators. 
One can check, using (2.25) and (2.1?:), that the mean 

and the variance of : 



W / x « 



are equal to x£L x ( Xn ) 
and cr, respectively. Therefore, the coordinate vectors 
\ x N-i> ■ ■ ■ ' ^o f > obtained through equation (2.25) for 
i = 1, 2, . . . A", a rc distributed according to the probabil- 



ity density ( 2.22 ). We not e in passing that the sequence 
of positions x n given by ( [2.25 ) defines a stochastic pro- 
cess, namely the well known Brownian diffusion process. 

Repeating J\f times the sampling of V the wave func- 
t ion $ f(x,t) can be determined according to ( [2.21 ) and 
( |2.14| ). Unfortunately, the algorithm outlined is imprac- 
tical since it provides only a route to calculate ^>(x,t) for 
a chosen time r, but no systematic method to obtain the 
ground state energy and wave function, which requires a 
description for r — * 00. 

Fortunately, the above so-called Monte Carlo algo- 
rithm can be improved and used to determine simulta- 
neously both Eq and </>o- The basic idea is to consider 
the wave function itself a probability density. This im- 
plies that the wave function should be a positive definite 
function, a constraint which limits the applicability of the 
suggested method. By sampling the initial wave function 
^(x, 0) at A/o points one generates as many Gaussian ran- 
dom walks which evolve in time according to ( p. 22 ) or, 
equivalently, according to (2.25); instead of tracing the 



motion of each random walk separately, one rather fol- 
lows the motion of a whole ensemble of random walks 
simultaneously. The advantage of this procedure is that 
one can sample the wave function of the system, through 
the actual position of the random walks and the products 
of the weights W along the corresponding trajectories, 
after each time step At. This procedure, as we explain 
below, also provides the possibility to readjust the value 
of En after each time step and to follow the time evolu- 
tion of the system for as many time steps as are needed 
to converge to the ground state wave function and energy 

The procedure, the so-called DMC method, interprets 
the integrand in (2.14), i.e., 



W(x N ) P{x n ,xn-i) ■■■ W(x 2 ) P(x 2 ,xi) 



process 27V process 2N—1 



process 4 process 3 



x W(xi) P(xi,x ) *(x ,0) , 

process 2 process 1 process 



(2.26) 



as a product of probabilities and weights to be mod- 
eled by a series of sequential stochastic processes 
0, 1, 2, . . . 2JV. We will explain now how these processes 
are described numerically. 

Initial State. The 0-th process describes "particles" 
(random walks) distributed according to the initial wave 
function \&(xo,0), which is typically chosen as a Dirac 
5— function 



#(z ,0) = S(x-x ) 



(2.27) 



where Xq is located in an area where the ground state 
of the quantum s ystem is expected to be large. The ini- 
tial distribution ( [2.27 ) is obtained by simply placing all 
particles initially at position xq. 

Diffusive Displacement. As we explained earlier, the 
successiv e posi tions x„_i, x n in (2.26) can be generated 
through (2.25). The Monte Carlo algorithm produces 
then the positions x\ — xq + crpi, x-i = Xi + opi-, 
etc by generating the series of random numbers p n , n = 
1,2,.... 

Birth— Death Process. Instead of accumulating the 
product of the weight factors W for each particle, it is 
more efficient numerically to replicate (see below) the 
particles after each time step with a probability propor- 
tional to W(x n ). In this way, after each time step At, 
the (unnormalized) wave function is given by a histogram 
of the spatial distribution of the particles. The calcula- 
tion of the wave function ^S(x,t) can be regarded as a 
simulated diffusion-reaction process of some imaginary 
particles. 

In the replication process each particle is replaced by 
a number of 



particles, where int[x] denotes the integer part of x and 
where u represents a random number uniformly dis- 
tributed in the interval [0, 1]. In case m n — the particle 
is deleted and one stops the diffusion process; this is re- 
ferred to as a "death" of a particle. In case of m n = 1 
the particle is unaffected and one continues with the next 
diffusion step. In case m n = 2, 3 one continues with the 
next diffusion step, but begins also a new series (in case 
of m n = 3 two new series) of diffusive displacements 
starting at the present location x n . This latter case is 
referred to as the " birth" of a particle (of two particles 
for m„ = 3). From ( 2.28J ) one can see that at most two 
new particles can be generated whereas one would ex- 
pect that for int[W / (x„) + u] > 4 three or more new 
series would be started. This limitation on the birth rate 
of the particles is necessary in order to avoid numeri- 
cal instabilities, especially at the beginning of the Monte 
Carlo simulation, when En may differ significantly from 
Eq. The error resulting from the limitation m n < 3 is 
expected to be small since for sufficiently small At holds 

W (x) « I - ^"^ Ar (2.29) 

n 

which, evidently, assumes values around unity. 



C. The Diffusion Monte Carlo Method 

We want to summarize now the algorithmic steps actu- 
ally taken in a straightforward application of the DMC 
method. To realize the suggested algorithm one starts 
with A/o "particles" , referred to as "replicas" , which are 
placed according to a distribution \I> (xo,0). As pointed 
out above, the actual numbers of replicas will vary as 
replicas "die" and new ones are "born". The replicas 
are characterized through a position x„ where the suf- 
fix n counts the diffusive displacements and where (j) 
counts the replicas. The number of replicas after n dif- 
fusive displacements will be denoted by J\f„. In the ini- 
tial step one samples No replicas assigning positions x , 
(i = 1, . ..,Nq) according to the distribution \& (xq,0). 
As mentioned, we actually choose all replicas to begin 
at the same point Xq. 

Rather than following the fate of each replica and of its 
descendants through many diffusive displacements, one 
follows all replicas simultaneously. Accordingly, one de- 
termines the positions xf following equation ( 2.25| ), i.e., 
one sets 



15) 



„U) 



v^aT/ 



GO 
m ■ p{ 



(2.30) 



for replicas j — 1,2,.. .No generating appropriate ran- 
dom numbers p\ . This can be regarded as a one-step 
diffusion proce ss of the system of replicas. Once the 
new positions ( |2.3C| ) have been determined, one evalu- 
ates the weight W(x^') given by ( 2.16 ) and, according 



[int[W(aj„) + u], 3] 



(2.28) to (2.28), one determines the set of integers 



,(J) 



for 



1,2, ...A/q. Replicas j with mf 



U) 



are ter- 



s.n 



minated, replicas ml — 1,2,3 are left unaffected, ex- 

cept, that in case m" = 2,3 one or two new replicas 
j' are added to the system and their position is set to 

Xi — xf . The number of replicas is counted and, 
thus, A/"i is determined. During the combined diffu- 
sion and birth-death process the distribution of replicas 
changes in such a way that the corresponding coordinates 
x\ J> will now be distributed according to a probability 
density which will be identified with ^>(x, At). 
Adaptation of the Reference Energy Er. As a re- 
sult of the birth-death process the total number Ai of 
replicas changes from its original A/o value. As discussed 
above, for Er values less then the ground state energy 
the distribution ^>(x,t) decays asymptotically to zero, 
i.e., the replicas eventually all die; for En values larger 
than the ground state energy, the distribution will in- 
crease indefinitely, i.e., the number J\f n will exceed all 
bounds. Only for Er = Eq can one expect a stable 
asymptotic distribution such that the number of replicas 
fluctuates around an average value Ao- The spatial distri- 
bution and increases/decreases of the number of replicas 
allow one to adjust the value of En as to keep the total 
number of replicas appro ximately constant. To this end, 
one proceeds from (2.29). Averaged over all replicas this 
equation reads 



(W)i 



1 



(Vh 



E 



R 



At 



where the average potential energy is 



„0")> 



A/'i 



(2.31) 



(2.32) 



i=i 



The reader should note that this average varies with the 
number of diffusion steps taken for all replicas; in the 
present case we consider the average after the first diffu- 
sive displacement. One would like the value of (W) to be 
eventually always unity such that the number of replicas 
remains constant. This leads to the proper choice of Er 



E 



(i) 



(Vh 



(2.33) 



However, if the distribution of replicas deviates too 
strongly from (f>o(x) the number of replicas can experi- 
ence strong changes which need to be repaired by 'over- 
compensating' through a suitable choice of Er. In fact, 
in case A/'i /A/o < 1 one would like to increase subse- 
quently the number of replicas in order to restore the 
initial number of replicas and, hence, choose an Er value 
larger than (2.33); in case Ai/A/o > 1 one would like to 



decrease subsequently the number of replicas and, hence, 



choose an Er value smaller than (2.33). A suitable choice 



E 



(2) 



^i + a7" 



M 

A/ 



(2.34) 



In order to justify that ( [2.34 ) serves our purpose one 
notes that exp[— At((V) — Er)/T%], for sufficiently small 
At, is 1 — At((V) — Er)/Ti. In the subsequent birth- 
death process this would lead to an A/i value (recall that 
during the diffusion process the number of replicas does 
not change) related to Er by 



*«-<*>-£0-& 



(2.35) 



However, one would like the total number of replicas dur- 
ing the next birth-death process to return to the initial 
val ue A/ p and Er to be given by an expression similar 
to ( 2.33| ). T his c an be achieved by adding to both sides 
of equation ( [2.35 ) the same quantity K (1 — Ai/Ao) /At. 
Hence, the redefined value for the reference energy is 






At 



A^ 
A/ 



(2.36) 



which 

(Pi 



by taking ( J2.33 ) into account, is identical to 



Equation (2.34) should be regarded as an empirical re- 
sult rather than an exact one. In fact, any expression 
of the form E R = (V) + a (1 - N/Nq), with arbitra: 



positive a, can be used equally well in place of (2.34)1 
Usually, the actual value of the "feedback" parameter a 
is chosen empirically for each individual problem so as 
to reduce as much as possible the statistical fluctuations 
in A/o and, at the same time, to diminish unwanted cor- 
relations b etwe en the successive generation of replicas. 
Equation ( [2.34 ) suggests that a good starting value for 
a is Ti/ At, a value used in our DMC program. 
2nd, 3rd, . . . Displacements. The diffusive displace- 
ments, birth-death processes and estimates for new Er 
values are repeated until, after a sufficiently large number 
of steps, the energy Er and the distribution of replicas 
becomes stationary. The ground state energy is then 



E: 



= lim (V) n 

n — >oo 



(2.37) 



since N n /N n -i — ► 1. The distribution of replicas pro- 
vides the ground state function 4>o(x). 

D. Systems with Several Degrees of Freedom 

The DMC method described above is valid for a quan- 
tum system with one degree of freedom x. However, the 
method can be easily extended to quantum systems with 
several degrees of freedom. Such systems arise, for exam- 
ple, in case of a particle moving in two or three spatial 
dimensions or in case of several interacting particles. For 
these two cases the DMC method can be readily gener- 
alized. 

Particle in d Dimensions. The Hamiltonian of a par- 
ticle of mass m moving in a potential V (x±, . . . Xd) can 
be written 



H 



2m *-^ dx 2 

a—l " 



V( Xl 



,x d ) 



(2.38) 



where x a , ct = 1, ..., d, denotes the Cartesian coordinates 
of the particle and d (= 2, 3) re presen ts the spatial di- 
mension. For the Hamiltonian ( 2. 38] ), the P MC algo- 
rithm can be devised in a similar way as for (2.2). The 
only difference to the one degree of freedom case is that 
during each time step At one needs to execute d ran- 
dom wal ks for each replica. Indeed, the kinetic energy 
term in ( 2.38| ) can be formally regarded as the sum of 
kinetic energies of d "particles" having the same mass m 
and moving along the x a , a = 1, ...,d, directions. Conse- 
quently, the diffusive displacements are governed by the 
distribution 



-* V*£n,l > X-n — 1,1 ! ■ ■ ■ ^n,d: Xfi— X,d) 

nLi(K7H?)*exp 



(2.39) 
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The product of probabilities can be described through 
independent ran dom processes, i.e., one can reproduce 
the probability ( 2.39| ) through d independent diffusive 
displacements applied for each replica to the d spatial 
directions. 

S Particles. In the case of S interacting particles, which 
move in d spatial dimensions, the most general form of 
the Hamiltonian is given by (assuming that no internal, 
e.g., spin, degrees of freedom are involved) 



H 



E 

i=i 



d 



ff 2 



2m, ^ dx 2 ja 



+ V({x ja }) 



(2.40) 



where V ({xj a }) accounts for both an interaction be- 
tween particles and for an interaction due to an external 
field; {xj a } denotes the dependence on the coo rdina tes 



of all particles. By rescaling the coordinates in ( 2.40J ) 



u ja 



jot ' 



j = 1,...,S , a = l,...,d 

(2.41) 



where m is a n arbitrary mass, one can make the Hamil- 
tonian (2.40) look formally the same as the Hamiltonian 
(2.38) describing a single particle of mass m which moves 
in d! = S x d spatial dimensions. Hence, the generaliza- 
tion of the DMC algorithm for this cas e is ag ain apparent. 
Sign Problem. The case in which ( 2. 40] ) actually de- 
scribes a system of identical particles cannot be treated 
like that of a single particle moving in d! = S x d spa- 
tial dimensions. For such systems a prescribed boson or 
fermion symmetry of the wave function with respect to 
an interchange of particles must be obeyed. For bosons 
(particles with integer spin) the total wave function (i.e., 
the product of the orbital and the spin wave functions) is 
symmetric with respect to any permutation of the parti- 
cles while for fermions (particles with half-integer spin) 



the total wave function is antisymmetric with respect to 
such permutation. This constraint determines the sym- 
metry of the orbital part of the ground state wave func- 
tion for fermionic systems (but not for bosonic ones): in 
many cases of interest, the (orbital) ground state wave 
function of fermionic systems will have nodes, i.e., re- 
gions with different signs, which make the DMC method, 
as presented here, inapplicable. We shall not address 
this issue in further detail and shall consider only cases 
where the so-called "sign problem" of the ground state 
wave function does not arise. 

To summarize, a DMC algorithm for a system with 
one degree of freedom can be adopted to a system of S 
interacting particles moving in d spatial dimensions with 
an effective dimension of d! = S x d. Exactly this feature 
of the DMC method makes it so attractive for the evalu- 
ation of the ground state of a quantum system. However, 
in the case of identical fermions, one needs to obey the 
actual symmetry of the ground state wave function and 
the method often is not applicable. 



III. ALGORITHM 

Our goal is to provide an algorithm for the DMC 
method presented in the previous section and to apply 
this algorithm to obtain the ground state energy and 
wave function for sample quantum systems. Some of 
the examples chosen below, e.g., the harmonic oscillator, 
have an analytical solution and, therefore, allow one to 
test the diffusion Monte Carlo method. Other examples, 
e.g., the hydrogen molecule, can not be solved analyti- 
cally and, hence, the DMC method provides a convenient 
way of solving the problem. The obtained results turn 
out to be in good agreement with_.results obtained by 
means of other numerical methodalj. 



A. Dimensionless Units 

In order to implement the DMC method into a numeri- 
cal algorithm, one needs to rewrite all the relevant equa- 
tions in dimensionless units. One can go from conven- 
tional (e.g., SI) units to dimensionless units by explicitly 
writing each physical quantity as its magnitude times the 
corresponding unit. In mechanics the unit of any phys- 
ical quantity can be expressed as a proper combination 
of three independent units, such as, L, T and £, which 
denote the unit of length, time and energy, respectively. 
By denoting the value of a given physical quantity with 
the same symbol as the quantity itself (e.g., xL — » x in 



the case of coordinate x), the Schrodinger equation (2.1C) 
can be recast 



~dr~ 



TiT d 2 f T£ 



V(x) - E R ] * 



2mL 2 dx 2 h 
It is convenient to choose L, T, and £ such that 



(3.1) 



fiT 
2mL 2 



and 



TS 

IT 



1 



(3.2) 



holds. Since there are three unknown units and only two 
relationships between them, one has the freedom to spec- 
ify the actual value of either L, T, £ while the value of 
the other two units follows from (3.2). 



In dimensionless units obeying (3.2) the original 
imaginary-time Schrodinger equation reads 






1<9 2 * 



2 dx 2 



[V{x) - E R ] * 



(3.3) 



It is not difficult to transcribe also the other relevant 
equations above in dimensionless units; for example, the 
functions (|2.15, 2.16) become 
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and 



W(x n ) = exp{-[V(x n )-E R }} 



respect ively. As a consequence, the diffusive displace- 
ments (2.25) are described by 



X n -l + V At/9„ 



(3.4) 



M max = 2000, r = 1000, At = 0.1, x„ 
Xmax = 20 and m, = 200. 



-20, 




Fig. 1. Flow diagram of the DMC algorithm. 



B. Computer Program 

The flow diagram of the computer programllj imple- 
menting the DMC method is shown in Fig. [y. Each 
block in the diagram performs specific tasks which are 
explained now. 

All external data required for the calculation are col- 
lected through a menu driven, interactive interface in the 
Input block. First, one has to select the quantum system 
on which the calculation is performed. At the program 
level this means to define the right spatial dimensionality 
d and the potential energy V (see Sec. [V) which corre- 
sponds to the selected system. The quantum systems 
covered by our program are the ground state of the har- 
monic oscillator, of the Morse oscillator, of the hydrogen 
atom, of the H \ ion (electronic state) and of the hydrogen 
molecule (electronic state) . The results of the simulation 
for all these cases are presented in Sec. IV below. The 
other input parameters are: the initial number of repli- 
cas (A/o), the maximum number of replicas (Af m ax), the 
seed value for the random number generators, the num- 
ber of time steps to run the simulation (tq), the value 
of the time step (At), the limits of the coordinates for 
the spatial sampling of the replicas [x m %ni x max ) and, 
finally, the number of spatial "boxes" (n&) for sorting 
the replicas during their sampling. Suggested values for 
these parameters are (in dimensionless units): A/q = 500, 



In the next block, Initialize replicas, a two di- 
mensional matrix, called psips, is defined with the fol- 
lowing structure: The first (row) index identifies the 
replicas and takes integer values between one and M max ■ 
The second (column) index points to information regard- 
ing the replica identified by the first index, and is a posi- 
tive integer less or equal to the num ber o f degrees of free- 
dom of the system (i.e., d! , see Sec. IID| ) plus one. For a 
given replica, say i, psips [i] [1] is used as an existence 
flag: it is zero (one) if the replica is dead (alive). The 
other elements psips [i] [j] , (J = 2, . . . , d' + 1) are used 
to store the coordinates of th e replica (i.e., Xj. a according 
to the notations of Sec. |ll D| ) during the simulation. 

To initialize the matrix psips one sets equal to one the 
value of the existence flag psips [j] [1] for j = 1, . . . A/o 
( replicas j > M a are not born yet and, accordingly, 
their existence flag is set to zero), and then one assigns 
the same coordinates for all these replicas (for the actual 
values of these coordinates, see Sec.pV|). Such a choice 
for the initial distribution of the replicas corresponds to 
a (5-function for ^(x, 0). For a suitable choice of xq one 
can be certain that there is always a significant overlap of 
vE^a;, 0) and the ground state wave function 4>o(x). The 
initial value of the reference energy E R is simply given 
by the value of the p oten tial energy at the initial position 
of the replicas [c.f. ( BJ33 )], 

After initializing the replicas the program enters into 



a loop which essentially consists of three routines: walk, 
branch and count. One loop corresponds to taking one 
time step At. 

walk: This routine performs the diffusion process of the 
replicas by adding to the coordinates of the active (alive) 
replicas the value vArp, where p is a random number 
drawn from a Gaussian distri buti on with mean zero and 
standard deviation one [c.f. (3.4)]. The program uses 
for thiSrfiurpose the Gaussian random number generator 
gasdevu . 

branch: The birth death (branching) processes, which 
follow the diffusion steps of the replicas, are performed 
by the routine branc h. For each alive replica the number 
m n , given by ( [2.28 ), is calculated. For the generation 
of u in ( [2.2SJ ), a uniformly distributed random number 
in the interval [0,1], the program employs the function 
ran3u. If a value m„ = results, the replica is killed by 
setting the corresponding existence flag to zero. If a value 
rn n = 1 results, the replica is left as it is. If a value n = 2 
results, the replica is duplicated (the first inactive replica 
in psips is set active with the same coordinates as the 
original replica). If a value n = 3 results, two identical 
copies of the replica are generated in psips. The reader 
may note that never more than two copies are born; this 
limitation is necessary to prevent the uncontrolled growth 
of the number of replicas which psips, due to its finite 
size, might not be able to accommodate. Such growth 
might occur when all the replicas are located in (almost) 
the same place and when one can expect that large val- 
ues m n can arise. Since (by choice) the replicas initially 
are located in one and the same point of the configura- 
tion space the first diffusion process does not spread the 
replicas far enough and, for certain initial positions, m n 
could then become large for all replicas. The algorithm 
also terminates if all the m n 's assume values zero. To 
avoid this possibility one needs to choose the initial lo- 
cation of the replicas with care. In general, any point 
where the ground state wave function is large is a good 
choice. 

count: The role of this routine is to return the ground 
state wave function of the system (i.e., the spatial dis- 
tribution of the replicas) at the end of the simulation. 
To this end, the spatial interval (x m i n ,x max ) is divided 
equally into rib "boxes" (sub-intervals) for each degree 
of freedom and then, by employing standard numerical 
methods^, one counts the distribution of replicas among 
these "boxes". The counting process starts after tq (in 
units of At) time steps when the system has already 
reached its stationary state (identified through a con- 
verged (V)t) and is performed in a cumulative way for 
another To time steps. This strategy can be justified 
as follows: once stationarity is reached, the wave func- 
tion \&(x, t) oc <fio( x ) will practically not change in time; 
hence, the replicas will essentially sample one and the 
same wave function at any subsequent time and the cu- 
mulative counting of the replicas in the "boxes" can be 
used; this procedure yields better statistics for (j>o(x) than 



sampling of replicas at only one instant in time; by cu- 
mulative counting, the effective number of replicas used 
to sample the wave function is enhanced by a factor of 
To (number of time steps the counting is done) . 

Once the spatial distribution of replicas is known, one 
can normalize the distribution to obtain the ground state 
wave function using 



4>o(%i) 
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(3.5) 



As should be apparent by now, until the completion of 
the algorithm the routines walk and branch are called 
2to times while the count routine is called only To times, 
namely, during the second half of the calculation. This 
is controlled by the block test shown in Fig. H. 

Finally, the Output block returns the results of the 
simulation. These results are: (i) the average value of 
the reference energy (Er) ~ Eq, which is calculated dur- 
ing the second part of the simulation (i.e., from To to 
2to) when the system is already stabilized; (ii) the corre- 
sponding standard deviation 5Er; (iii) the (imaginary) 
time evolution of (En) for the first t time steps (used 
basically to check how fast stationarity is reached by the 
system during the simulation); (iv) the normalized spa- 
tial distribution of the replicas, i.e., the ground state wave 
function. Note that the average reference energy after n 
time steps is defined through 



1 n 
(Er{t = tiAt)) = -YE R (iAr) 



(3.6) 



IV. EXAMPLES 

In this section we report on the results obtained, by 
means of the DMC program, for the ground state energy 
and wave function of some quantum mechanical systems. 
The program was executed on an HP-9000 (series 700) 
workstation. In each case we specify the units used and 
present numerical results of the simulation. For all sim- 
ulat ions t he values of the input parameters suggested in 
Sec. IIIB, have been employed. 



A. Harmonic Oscillator 

For a one-dimensional harmonic oscillator, character- 
ized by the proper angular frequency u>, one has (in SI 

units) 



V(x) 



1 



2 *? 

rauj x 



(4.1) 
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Choosing T = 1/w, one obtains from ( |3.2| ) L — y'ti/mu 
and £ = Tilo. In corresponding dimensionless units, the 
exact ground state eaergy is E — 1/2 and the ground 
state wave function ia 



<Po(x) 



exp 



x 

~2 



(4.2) 



The results of the Monte Carlo simulation for the har- 
monic oscillator are given in Fig. g. At the beginning of 
the simulation all replicas are located at the origin. The 
main graph shows the rapid convergence of the reference 
energy (Er(t)) towards the exact ground state energy 
Eq . The inset contains the plot of the ground state wave 
function: the results of the simulation are represented by 
triangles, while t he c ontinuous line corresponds to the 
analytical result ( |4.2| ). The agreement between the dif- 
fusion Monte Carlo result and the analytical expressions 
is very good. 
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Fig. 2. Reference energy Er (in dimensionless units) as a 
function of the imaginary time r (in units of time steps At) 
obtained by the DMC method for the harmonic oscillator. 
Er converges rapidly towards the exact ground state energy 
Eo — 0.5. The inset shows the corresponding ground state 
wave function <j>o(x). The result of the simulation is repre- 
sented by dots while the continuous line corresponds to the 
analytical solution (4.2) 



B. Morse Oscillator 



The Morse potential is defined through 



V(x) = D(e 



-lax 



2e~ 



(4.3) 



In this case one has a natural length scale through 
L = 1/a whi ch c an be used as the unit of length. As a 
result, from (3.2), one has T — m/ha 2 and £ = h a 2 /m. 
For simplicity we shall consider only the case when (in 
dimensionless units) D = 1/2; the exact ground state en- 
ergy is E = —1/8, and the corresponding wave function 



M x ) = \/2cxp(-e x - ~) 



(4.4) 



The results of the DMC description are presented in 
Fig. 0. Initially, all replicas were positioned in our simu- 
lation at the origin. The figure represents the time evo- 
lution of (En) towards the exact ground state energy 



Eo = —0.125. The figure demonstrates also that the 
the resulting ground state wave function is in excellent 
agreement with (4.4). 
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Fig. 3. DMC description of the Morse oscillator. The time 
evolution of Er is shown. The inset presents the ground state 
wave function; the numerical result is represe nted through 
dots, the wave function in the analytical form (4.4) is repre- 
sented by a continuous line. 



C. Hydrogen Atom 

In case of the hydrogen atom it is customary to choose 
the unit of length equal to the Bohr radius a = h /me 2 (~ 



0.53A). Thus, by setting L — a, one finds T = h /me 
and £ = me 4 /h (ss 27.21eV r ). The well-known ground 
state energy (in dimensionless units) is Ep — —1/2 and 
the corresponding radial wave function istl 



Mr) = 2e- 



(4.5) 



We have carried out a DMC simulation for the hydro- 
gen atom, generalizing the previous description to three 
spatial dimensions. At r = all the replicas were located 
at (0, 0, 1), i.e., at a distance of one Bohr radius from the 
origin, along the positive z axis. Figure [| shows the con- 
vergence of (Er) to the exact ground state energy. This 
convergence, however, is not as rapid as in the case of the 
harmonic oscillator and the Morse oscillator. This is not 
surprising since the hydrogen atom has three degrees of 
freedom which require more sampling. The running time 
of the simulation has also increased (see Table ||) . The in- 
set shows both the radial wave function (j>o(r) (triangles) 
and the function x( r ) — r ' ( / l o( r ) (squares). For compari- 
son, the corresponding analytical solutions arc also plot- 
ted with continuous lines. Again the agreement between 
the analytical results and those obtained by our algo- 
rithm are good. The error of the radial wave function in 
the vicinity of the origin is due to insufficient sampling of 
the number of replicas in this region of the configuration 
space. To improve the wave function for r ~ one may 
decrease the size of the counting "boxes" in the vicinity 
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of the origin, increase the number of replicas, or increase 
r. 
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Fig. 4. DMC description of the hydrogen atom. The time 
evolution of En is shown. The inset presents both the radial 
wave function 4>o(r) and the function x( r ) — r( t > o( r )- 



D. HJ Ion 



The H^" ion is held stabilized through a single electron 
moving in the electric field of two protons separated, at 
equilibrium, by a distance R — 1.06^4tj. The ground 
state of this quantum system can not be determined an- 
alytically. The numerically determined ground state en- 
ergy of H+ is E = -16.252 ± 0.002 eVt3. The DMC 
method can be applied in this case as for the hydrogen 
atom. The Hamiltonian is 



H 2 „ 

H = V 2 

2m 



\R\ \r - \R\ 



(4.6) 



where R denotes the separation between the two pro- 
tons. The results of the DMC description are presented 
in Fig. pi The same dimensionless units as in the case 
of the hydrogen atom were employed and replicas were 
located initially at the origin, the nuclei being located 
at (0,0, ±12/2), for R = 2. Figure | shows that the 
ground state energy obtained asymptotically is —16.75 
eV (see Table_|) which differs from the more exact nu- 
merical valueEJ by 0.5 eV, i.e., by about 3%. The inset 
in Fig. H shows a plot of the spatial distribution of the 
replicas, demonstrating that the electronic ground state 
wave function is nearly spherically symmetric. 



E. H 2 Molecule 



The H2 molecule is formed by 
librium separated by R — 0.74AI 



o protons, at equi- 
and by two elec- 



trons. In the Born approximation one considers the pro- 
ton positions fixed and solves the corresponding station- 
ary Schrodinger equation for the electrons. The wave 
function of the two electrons must be antisymmetric with 
respect to exchange of the electrons. In the ground state 
the electrons are in a singlet spin state 
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Fig. 5. DMC description of the Hj ion. The time evolution 
of Er is shown. The inset presents the spatial distribution of 
the replicas, i.e., the electron cloud. 



x(l>2) = y-[xi, + |(l)xi,-i(2)-Xi > -i(l)Xi, + i(2) 

(4.7) 

which is antisymmetric. Here xi ± i(1,2) denotes the 
wave function of electron 1,2 in the spin i state with 
magnetic quantum numbers ±^. Accordingly, in the 
wave function of the electronic ground state 



*o = fc(ri,r 2 ) X (l,2) 



(4.8) 



the factor <&(ri,r2), describing the spatial distribution 
of the electrons, must be symmetric. This factor can be 
determined through 
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where 4>{r\,r2) is a solution of 
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TABLE I. Results obtained during two different simulations performed on four quantum systems listed in the first column. 
During Simulation I (II) the initial number of replicas was 500 (4000), the time step At = 0.1 (0.05), the length of simulation 
r = 1000 Ar (2000 Ar), the random number generator seed 1 and the number of boxes used to calculate the spatial distribution 
of the replicas 200 (400). {En) and SEr are defined in the text. At represents the actual running time of the simulation on 
a HP-9000 (series 700) workstation. For the hydrogen atom the energies are given both in dimensionless units and in eV (in 
parenthesis), respectively. For comparison, Eq in the last column represents the exact value (analytical or obtained through an 
alternative numerical method) of the corresponding ground state energy. 







Simulation I 






Simulation II 






Quantum system 


{En) 


5E R 


At 


{Er) 


SE R 


At 


Eq 


Harmonic oscillator 


0.505 


0.094 


8 sec 


0.500 


0.048 


4 min 


0.5 


Morse oscillator 


-0.1236 


0.0749 


10 sec 


-0.1245 


0.0330 


4 min 


-0.125 


Hydrogen atom 


-0.495 


0.080 


18 sec 


-0.505 


0.040 


5 min 


0.5 




(-13.477 eV) 


(2.186 eV) 




(-13.752 eV) 


(1.093 eV) 




(13.6 eV) 


Hj ion 


-16.753 eV 


2.869 eV 


40 sec 


-16.476 eV 


1.389 eV 


11 min 


-16.25(2) eV 


H2 molecule 


-30.973 eV 


3.638 eV 


55 sec 


-31.968 eV 


1.754 eV 


16 min 


-31.6(87) eV a 



a Numerical values calculated based upon the heat of dissociation given by Hertzber 
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The Schrodinger equation (4.10) cannot be solved an- 
alytically. The available numerical result for the-ground 
state energy of H 2 is E = -31.688 ± 0.013 e\M The 
results of the diffusion Monte Carlo simulation for H2 are 
shown in Fig. 0. The same dimensionless units as in the 
case of the hydrogen atom were employed and the repli- 
cas at r = were located initially at {(0, 0, 1)(0, 0, —1)}; 
the position of the protons were (0,0,±i?/2), with R = 
1.398. The simulation took about 55 seconds to com- 
plete and the obtained ground state energy was E = 
-30.973 eV (see Table ||. This result differs from the ex- 
act energy by 0.715 eV, i.e., by 2.3%. The inset of Fig. | 
shows the spatial distribution of the replicas which is 
nearly spherically symmetric. In Fig. p] we present the 
results of a calculation in which an unphysically large R 
value of 8.398 was assumed. In this case the distribu- 
tion of the electronic cloud around the protons is clearly 
anisotropic, the energy of the electrons is still negative. 




Fig. 7. En vs. r and the spatial distribution of the replicas 
(inset) at the end of the simulation for the H2 molecule for 
an unphysically large separation between the two protons of 
8.398 Bohr radius. 
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Fig. 6. DMC description of the electronic ground state of 
the H2 molecule. The time evolution of Er is shown. The 
inset presents the spatial distribution of the replicas, i.e., the 
electron cloud. 



V. DISCUSSION 

In this article, we have presented a detailed account 
of the path integral formulation of the DMC method 
devised for calculating simultaneously the ground state 
energy and wave function of an arbitrary quantum sys- 
tem. A simple numerical algorithm based on the DMC 
method has been formulated and a computer program, 
based on this algorithm, has been applied to determine 
numerically the ground state of a few quantum systems 
of pedagogical interest. 

The DMC algorithm, as presented in this article, is 
quite unstable numerically. We want to demonstrate here 
the need for an improvement of the method. This need 
arises due to the strong fluctuations stemming from the 
birth-death processes employed in this method. These 
fluctuations affect the reference energy Er(t) which is 
expected to converge to the ground state energy. This 
convergence is observed only for the value of {Er) defined 
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through equation (3.6). The actual time dependence of 
Er is presented in Fig. g for the case of the harmonic 
oscillator. The fluctuations in Er are large, but symmet- 
rically distributed around (En) which explains the good 
agreement between the result of the simulation with the 
exact result. Attempts to use our DMC program to cal- 
culate the electronic ground state energy of the Hj ion 
and H2 molecule as a function of the separation between 
the two protons, which would allow one to determine the 
equilibrium separation from the minimum of this depen- 
dence, have failed due to large fluctuations in Er. 

The DMC algorithm can be significantly improved^! 
sorting to a method called "importance sampling" Erl 
The basic idea of this method is to change the prob- 
ability distribution of the replicas in a controlled way. 
This can be achieved by reformulating (2.10) such that 
the resulting equation has a solution ^(x, t) multiplied 
by an approximation of the ground state wave function, 
the latter being obtained, for example, from a variational 
method. Application of the DMC method to this new 
equation yields then replicas which spend more time in 
"important" regions of the configuration space where the 
wave function ^(x : r) is expected to be large. For details 
regarding importaufitfpSampling the reader is referred to 
the literature cite^ 1 ' 1 ' 



Any implementation of the DMC method leads to sys- 
tematic errors due to the use of a finite time step At. 
Apparently these errors can be reduced to zero by choos- 
ing At — > 0. Unfortunately, this is not the case. Even 
worse: making the time step shorter and shorter does not 
only increase the needed computer time, but also renders 
the successive generations of replicas more and more cor- 
related such that their distribution actually departs more 
and more from the ground state wave function. For each 
quantum system investigated it is necessary to find the 
most convenient value of At which is short enough to 
produce small systematic errors and at the same time is 
long enough to keep the successive distributions of repli- 
cas sufficiently uncorrelated. 



Fig. 8. Fluctuations of En about its average value 
(Er) ~ 0.5 well after the distribution of replicas reached its 
stationary form. The relatively large fluctuations are due to 
the birth-death processes which persist even if the replicas are 
distributed according to the exact ground state wave function 
</>o(x) and the reference energy En has been adapted to the 
exact ground state energy. 

In order to apply the DMC method to calculate the 
ground state properties of a system with interacting 
fcrmions one has to treat, as mentioned above, the "sign 
problem" due to the antisymmetry property of the many- 
fermion wav&fuiiction. Two principle methods, jthja fixed- 
node methoduJ'u and the release-node methodotl, have 
been proposed to deal with this problem. 

Once some kind of importance sampling is imple- 
mented and the sign-problem, in the case of many- 
fermion systems, is resolved, the DMC method can be 
used to compute ground state properties for molecules 
or molecular clusters^! and for quantum spin, boson and 
fermion systemaii The method is also applicable to the 
study of ground-state phase transitions due to quan- 
tum fluctuations, a topic of modern condensed matter 
physics^. 

Finally we would like to mention that the DMC 
method has been extended and successfully applied ta 
the study of the excited states of molecules and clusterstl 
and also to the study of finite temperature properties 
of different condensed matter systems. Also, the DMC 
method._b.as been successfully applied in quantum field 
theories!!!! 
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